library(MendelianRandomization)
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
library(ks)
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2     ✓ purrr   0.3.4
✓ tibble  3.0.3     ✓ stringr 1.4.0
✓ tidyr   1.1.0     ✓ forcats 0.5.0
✓ readr   1.3.1     
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(stringr)

First let’s load our data into a dataframe. This dataframe has 5 non-null columns: 1. seq_num: Index of the DNA sequence from which the effect sizes and standard errors come. 2. X_pred_mean: Effect size for the \(Z \rightarrow X\) relationship. Determined by taking the difference between predicted transcription factor binding probability for a reference sequence and a mutated version of the reference. 3. X_pred_var: Standard error of the \(Z \rightarrow X\) effect size estimate. 4. Y_pred_mean: Effect size for the \(Z \rightarrow Y\) relationship. Determined by taking the difference between predicted chromatin accessibility probability for a reference sequence and a mutated version of the reference. 5. Y_pred_var: Standard error of the \(Z \rightarrow Y\) effect size estimate.

data_dir <- "../../dat/"
results_dir <- "../../dat/"
exposure_tf <- params$exposure_tf
outcome_tf <- params$outcome_tf

results_fname <- str_c(exposure_tf, outcome_tf, "effect_sizes.csv", sep = "_")
seq_predictions <- read.csv(file.path(results_dir, results_fname))
is.nan.data.frame <- function(x) {
  do.call(cbind, lapply(x, is.nan))
}
seq_predictions = seq_predictions %>% drop_na()

n_seqs <- max(seq_predictions["seq_num"])
print(str_c("Running on ", n_seqs, " sequences' predictions for exposure ", exposure_tf, " and outcome ", outcome_tf))
[1] "Running on 60 sequences' predictions for exposure Nanog and outcome Sox2"
seq_predictions
ivw_vals <- list()
egger_vals <- list()

egger_result <- matrix(ncol = 8, nrow = n_seqs)
plots <- c()

for (seq in (1:n_seqs)) {
  seq_i_predictions <- subset(seq_predictions, seq_num == seq)

  bxt <- unlist(seq_i_predictions["X_pred_mean"])
  bxset <- unlist(seq_i_predictions["X_pred_var"])
  byt <- unlist(seq_i_predictions["Y_pred_mean"])
  byset <- unlist(seq_i_predictions["Y_pred_var"])
  MRInputObject <- mr_input(
    bx = bxt,
    bxse = bxset,
    by = byt,
    byse = byset
  )
  EggerObject <- mr_egger(MRInputObject)
  plots <- c(plots, mr_plot(MRInputObject))
  egger_result[seq, ] <- c(
    EggerObject$Estimate,
    EggerObject$CILower.Est,
    EggerObject$CIUpper.Est,
    EggerObject$I.sq,
    EggerObject$Pleio.pval,
    EggerObject$StdError.Est,
    EggerObject$Intercept,
    EggerObject$Pvalue.Est
  )
}

colnames(egger_result) <- c(
  "ce",
  "cil",
  "ciu",
  "i.sq",
  "pleio",
  "std",
  "int",
  "pval"
)
egger_result <- as_tibble(egger_result)
# egger_result$int
# egger_result$pleio
egger_result$pval
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
egger_result$cil
 [1] 0.6624111 0.4946017 0.2596899 0.3392746 0.2643942 0.3214277 0.3202355 0.4878342 0.3956594 0.4714063 0.4591280 0.4529287 0.1941315 0.3179280 0.4748364 0.2723454 0.3476054 0.3609901 0.5434091 0.3501268 0.3697978 0.1090658 0.3051474 0.3916034 0.3526319
[26] 0.3088198 0.5313684 0.2309187 0.6876855 0.5278499 0.2404396 0.3849298 0.4827181 0.5442420 0.6083380 0.5331485 0.3213684 0.3346143 0.3504754 0.4524046 0.3085015 0.5211029 0.5249166 0.3742332 0.4465752 0.4628427 0.5612626 0.4826918 0.5524170 0.5218692
[51] 0.5326102 0.4124616 0.3531925 0.4569031 0.3447336 0.2388891 0.5691459 0.4368713 0.4230229 0.5664108
egger_result$ciu
 [1] 0.6925329 0.5313000 0.2936519 0.3674658 0.2868699 0.3560291 0.3481173 0.5225789 0.4305670 0.5033834 0.4898804 0.4843012 0.2326947 0.3654092 0.5080045 0.3106776 0.3828946 0.4023319 0.5752519 0.3835175 0.4043459 0.1375472 0.3386546 0.4217517 0.3915041
[26] 0.3429232 0.5702170 0.2619052 0.7127931 0.5620440 0.2704281 0.4131669 0.5121622 0.5792549 0.6412172 0.5644271 0.3477196 0.3631897 0.3827492 0.4843635 0.3397468 0.5549002 0.5585180 0.4082484 0.4824702 0.4942072 0.5997734 0.5204230 0.5942819 0.5529227
[51] 0.5619741 0.4438768 0.3866925 0.4862023 0.3754247 0.2682767 0.6059248 0.4702945 0.4504525 0.5946959

Now we have the causal effect estimates, confidence intervals, and (for Egger) pleiotropy p-values.

egger_result

Let’s summarize the IVW & Egger causal effect estimates as histograms. This will give us an idea of how heterogeneous the alleged causal relationships are at different regions of the genome.

d <- density(egger_result$ce)
ce_smoothed <- ks::kde(
  x = egger_result$ce,
  binned = TRUE,
  compute.cont = TRUE,
  xmin = c(min(egger_result$ce) - 1),
  xmax = c(max(egger_result$ce) + 1),
  bgridsize = c(200)
)

if (params$save_plots) {
  output_fname <- str_c(cell_type, tf, "ces_kde.png", sep = "_")
  output_fpath <- file.path(params$output_dir, output_fname)
  png(output_fpath, width = 512, height = 320)
  par(mar = c(5, 6, 5, 1) + .1)
}
plot(
  ce_smoothed,
  xlab = "Causal Effect",
  # main = str_c("Causal Effect Estimates (", tf, ")"),
  cex.lab = 2,
  cex.main = 2,
  ylab = "Density"
)

if (params$save_plots) {
  dev.off()
}
sorted_idxs <- sort(egger_result$ce, index.return=TRUE)$ix
median_idx <- sorted_idxs[length(sorted_idxs) / 2 + 1]
seq_median_preds <- subset(seq_predictions, seq_num == median_idx)

bxt <- unlist(seq_median_preds["X_pred_mean"])
bxset <- unlist(seq_median_preds["X_pred_var"])
byt <- unlist(seq_median_preds["Y_pred_mean"])
byset <- unlist(seq_median_preds["Y_pred_var"])

MRInputObject <- mr_input(
  bx = bxt,
  bxse = bxset,
  by = byt,
  byse = byset
)
MedianEggerObj <- mr_egger(MRInputObject)
mr_plot(MRInputObject, line="egger")

if (params$save_plots) {
  png(file.path(params$output_dir, output_fname), width = 512, height = 320)
}
par(mfrow=c(2, 2))

for (seq in (1:n_seqs)) {
  seq_i_predictions <- subset(seq_predictions, seq_num == seq)

  bxset <- unlist(seq_i_predictions["X_pred_var"])
  byset <- unlist(seq_i_predictions["Y_pred_var"])
  
  plot(bxset, byset, main = str_c("CA vs. TF std error"))
}

if (params$save_plots) {
  dev.off()
}
dmode <- function(x) {
  den <- density(x, kernel = c("gaussian"))
  (den$x[den$y == max(den$y)])
}
dmode(egger_result$ce)
median(egger_result$ce)
length(egger_result$ce[egger_result$ce < 0])
length(egger_result$ce)

Let’s also look at the distribution of pleiotropy p-values for different sequences. If they’re concentrated around .5, this means that the algorithm believes there to be very little pleiotropy.

qplot(pleio, data = egger_result, binwidth = .04)
p <- ggplot2::ggplot(egger_result, aes(x = ce, y = std)) +
  geom_point()
p
if (params$save_plots) {
  output_fname <- str_c(cell_type, tf, "ex_eff_sizes.png", sep = "_")
  output_fpath <- file.path(params$output_dir, output_fname)
  png(output_fpath, width = 512, height = 320)
  par(mar = c(5, 6, 5, 1) + .1)
}

plot(bxt, byt, xlab=str_c(exposure_tf, " Effect Sizes"), ylab=str_c(outcome_tf, " Effect Sizes"))

if (params$save_plots) {
  dev.off()
}
LS0tCnRpdGxlOiAiTWVuZGVsaWFuIFJhbmRvbWl6YXRpb24gYW5hbHlzaXMgb2YgRGVlcFNFQSBnZW5lcmF0ZWQgZGF0YSIKb3V0cHV0OiBodG1sX25vdGVib29rCnBhcmFtczoKICBleHBvc3VyZV90ZjogIktsZjQiCiAgb3V0Y29tZV90ZjogIk9jdDQiCiAgb3V0cHV0X2RpcjogIi4uLy4uL2ZpZy9SLyIKICBzYXZlX3Bsb3RzOiBGCi0tLQpgYGB7cn0KbGlicmFyeShNZW5kZWxpYW5SYW5kb21pemF0aW9uKQpsaWJyYXJ5KGRwbHlyKQpsaWJyYXJ5KGtzKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KHN0cmluZ3IpCmBgYAoKRmlyc3QgbGV0J3MgbG9hZCBvdXIgZGF0YSBpbnRvIGEgZGF0YWZyYW1lLiBUaGlzIGRhdGFmcmFtZSBoYXMgNSBub24tbnVsbCBjb2x1bW5zOiAKMS4gYHNlcV9udW1gOiBJbmRleCBvZiB0aGUgRE5BIHNlcXVlbmNlIGZyb20gd2hpY2ggdGhlIGVmZmVjdCBzaXplcyBhbmQgc3RhbmRhcmQgZXJyb3JzIGNvbWUuCjIuIGBYX3ByZWRfbWVhbmA6IEVmZmVjdCBzaXplIGZvciB0aGUgXCggWiBccmlnaHRhcnJvdyBYIFwpIHJlbGF0aW9uc2hpcC4gRGV0ZXJtaW5lZCBieSB0YWtpbmcgdGhlIGRpZmZlcmVuY2UgYmV0d2VlbiBwcmVkaWN0ZWQgdHJhbnNjcmlwdGlvbiBmYWN0b3IgYmluZGluZyBwcm9iYWJpbGl0eSBmb3IgYSByZWZlcmVuY2Ugc2VxdWVuY2UgYW5kIGEgbXV0YXRlZCB2ZXJzaW9uIG9mIHRoZSByZWZlcmVuY2UuCjMuIGBYX3ByZWRfdmFyYDogU3RhbmRhcmQgZXJyb3Igb2YgdGhlIFwoIFogXHJpZ2h0YXJyb3cgWCBcKSBlZmZlY3Qgc2l6ZSBlc3RpbWF0ZS4KNC4gYFlfcHJlZF9tZWFuYDogRWZmZWN0IHNpemUgZm9yIHRoZSBcKCBaIFxyaWdodGFycm93IFkgXCkgcmVsYXRpb25zaGlwLiBEZXRlcm1pbmVkIGJ5IHRha2luZyB0aGUgZGlmZmVyZW5jZSBiZXR3ZWVuIHByZWRpY3RlZCBjaHJvbWF0aW4gYWNjZXNzaWJpbGl0eSBwcm9iYWJpbGl0eSBmb3IgYSByZWZlcmVuY2Ugc2VxdWVuY2UgYW5kIGEgbXV0YXRlZCB2ZXJzaW9uIG9mIHRoZSByZWZlcmVuY2UuCjUuIGBZX3ByZWRfdmFyYDogU3RhbmRhcmQgZXJyb3Igb2YgdGhlIFwoIFogXHJpZ2h0YXJyb3cgWSBcKSBlZmZlY3Qgc2l6ZSBlc3RpbWF0ZS4KCmBgYHtyfQpkYXRhX2RpciA8LSAiLi4vLi4vZGF0LyIKcmVzdWx0c19kaXIgPC0gIi4uLy4uL2RhdC8iCmV4cG9zdXJlX3RmIDwtIHBhcmFtcyRleHBvc3VyZV90ZgpvdXRjb21lX3RmIDwtIHBhcmFtcyRvdXRjb21lX3RmCgpyZXN1bHRzX2ZuYW1lIDwtIHN0cl9jKGV4cG9zdXJlX3RmLCBvdXRjb21lX3RmLCAiZWZmZWN0X3NpemVzLmNzdiIsIHNlcCA9ICJfIikKc2VxX3ByZWRpY3Rpb25zIDwtIHJlYWQuY3N2KGZpbGUucGF0aChyZXN1bHRzX2RpciwgcmVzdWx0c19mbmFtZSkpCmlzLm5hbi5kYXRhLmZyYW1lIDwtIGZ1bmN0aW9uKHgpIHsKICBkby5jYWxsKGNiaW5kLCBsYXBwbHkoeCwgaXMubmFuKSkKfQpzZXFfcHJlZGljdGlvbnMgPSBzZXFfcHJlZGljdGlvbnMgJT4lIGRyb3BfbmEoKQoKbl9zZXFzIDwtIG1heChzZXFfcHJlZGljdGlvbnNbInNlcV9udW0iXSkKcHJpbnQoc3RyX2MoIlJ1bm5pbmcgb24gIiwgbl9zZXFzLCAiIHNlcXVlbmNlcycgcHJlZGljdGlvbnMgZm9yIGV4cG9zdXJlICIsIGV4cG9zdXJlX3RmLCAiIGFuZCBvdXRjb21lICIsIG91dGNvbWVfdGYpKQpzZXFfcHJlZGljdGlvbnMKYGBgCgoKYGBge3J9Cml2d192YWxzIDwtIGxpc3QoKQplZ2dlcl92YWxzIDwtIGxpc3QoKQoKZWdnZXJfcmVzdWx0IDwtIG1hdHJpeChuY29sID0gOCwgbnJvdyA9IG5fc2VxcykKcGxvdHMgPC0gYygpCgpmb3IgKHNlcSBpbiAoMTpuX3NlcXMpKSB7CiAgc2VxX2lfcHJlZGljdGlvbnMgPC0gc3Vic2V0KHNlcV9wcmVkaWN0aW9ucywgc2VxX251bSA9PSBzZXEpCgogIGJ4dCA8LSB1bmxpc3Qoc2VxX2lfcHJlZGljdGlvbnNbIlhfcHJlZF9tZWFuIl0pCiAgYnhzZXQgPC0gdW5saXN0KHNlcV9pX3ByZWRpY3Rpb25zWyJYX3ByZWRfdmFyIl0pCiAgYnl0IDwtIHVubGlzdChzZXFfaV9wcmVkaWN0aW9uc1siWV9wcmVkX21lYW4iXSkKICBieXNldCA8LSB1bmxpc3Qoc2VxX2lfcHJlZGljdGlvbnNbIllfcHJlZF92YXIiXSkKICBNUklucHV0T2JqZWN0IDwtIG1yX2lucHV0KAogICAgYnggPSBieHQsCiAgICBieHNlID0gYnhzZXQsCiAgICBieSA9IGJ5dCwKICAgIGJ5c2UgPSBieXNldAogICkKICBFZ2dlck9iamVjdCA8LSBtcl9lZ2dlcihNUklucHV0T2JqZWN0KQogIHBsb3RzIDwtIGMocGxvdHMsIG1yX3Bsb3QoTVJJbnB1dE9iamVjdCkpCiAgZWdnZXJfcmVzdWx0W3NlcSwgXSA8LSBjKAogICAgRWdnZXJPYmplY3QkRXN0aW1hdGUsCiAgICBFZ2dlck9iamVjdCRDSUxvd2VyLkVzdCwKICAgIEVnZ2VyT2JqZWN0JENJVXBwZXIuRXN0LAogICAgRWdnZXJPYmplY3QkSS5zcSwKICAgIEVnZ2VyT2JqZWN0JFBsZWlvLnB2YWwsCiAgICBFZ2dlck9iamVjdCRTdGRFcnJvci5Fc3QsCiAgICBFZ2dlck9iamVjdCRJbnRlcmNlcHQsCiAgICBFZ2dlck9iamVjdCRQdmFsdWUuRXN0CiAgKQp9Cgpjb2xuYW1lcyhlZ2dlcl9yZXN1bHQpIDwtIGMoCiAgImNlIiwKICAiY2lsIiwKICAiY2l1IiwKICAiaS5zcSIsCiAgInBsZWlvIiwKICAic3RkIiwKICAiaW50IiwKICAicHZhbCIKKQplZ2dlcl9yZXN1bHQgPC0gYXNfdGliYmxlKGVnZ2VyX3Jlc3VsdCkKYGBgCgpgYGB7cn0KZWdnZXJfcmVzdWx0JGkuc3EKIyBlZ2dlcl9yZXN1bHQkcGxlaW8KZWdnZXJfcmVzdWx0JHB2YWwKZWdnZXJfcmVzdWx0JGNpbAplZ2dlcl9yZXN1bHQkY2l1CmBgYAoKTm93IHdlIGhhdmUgdGhlIGNhdXNhbCBlZmZlY3QgZXN0aW1hdGVzLCBjb25maWRlbmNlIGludGVydmFscywgYW5kIChmb3IgRWdnZXIpIHBsZWlvdHJvcHkgcC12YWx1ZXMuCmBgYHtyfQplZ2dlcl9yZXN1bHQKYGBgCgpMZXQncyBzdW1tYXJpemUgdGhlIElWVyAmIEVnZ2VyIGNhdXNhbCBlZmZlY3QgZXN0aW1hdGVzIGFzIGhpc3RvZ3JhbXMuIFRoaXMgd2lsbCBnaXZlIHVzIGFuIGlkZWEgb2YgaG93IGhldGVyb2dlbmVvdXMgdGhlIGFsbGVnZWQgY2F1c2FsIHJlbGF0aW9uc2hpcHMgYXJlIGF0IGRpZmZlcmVudCByZWdpb25zIG9mIHRoZSBnZW5vbWUuCmBgYHtyfQpkIDwtIGRlbnNpdHkoZWdnZXJfcmVzdWx0JGNlKQpjZV9zbW9vdGhlZCA8LSBrczo6a2RlKAogIHggPSBlZ2dlcl9yZXN1bHQkY2UsCiAgYmlubmVkID0gVFJVRSwKICBjb21wdXRlLmNvbnQgPSBUUlVFLAogIHhtaW4gPSBjKG1pbihlZ2dlcl9yZXN1bHQkY2UpIC0gMSksCiAgeG1heCA9IGMobWF4KGVnZ2VyX3Jlc3VsdCRjZSkgKyAxKSwKICBiZ3JpZHNpemUgPSBjKDIwMCkKKQoKaWYgKHBhcmFtcyRzYXZlX3Bsb3RzKSB7CiAgb3V0cHV0X2ZuYW1lIDwtIHN0cl9jKGNlbGxfdHlwZSwgdGYsICJjZXNfa2RlLnBuZyIsIHNlcCA9ICJfIikKICBvdXRwdXRfZnBhdGggPC0gZmlsZS5wYXRoKHBhcmFtcyRvdXRwdXRfZGlyLCBvdXRwdXRfZm5hbWUpCiAgcG5nKG91dHB1dF9mcGF0aCwgd2lkdGggPSA1MTIsIGhlaWdodCA9IDMyMCkKICBwYXIobWFyID0gYyg1LCA2LCA1LCAxKSArIC4xKQp9CnBsb3QoCiAgY2Vfc21vb3RoZWQsCiAgeGxhYiA9ICJDYXVzYWwgRWZmZWN0IiwKICAjIG1haW4gPSBzdHJfYygiQ2F1c2FsIEVmZmVjdCBFc3RpbWF0ZXMgKCIsIHRmLCAiKSIpLAogIGNleC5sYWIgPSAyLAogIGNleC5tYWluID0gMiwKICB5bGFiID0gIkRlbnNpdHkiCikKaWYgKHBhcmFtcyRzYXZlX3Bsb3RzKSB7CiAgZGV2Lm9mZigpCn0KYGBgCgpgYGB7cn0Kc29ydGVkX2lkeHMgPC0gc29ydChlZ2dlcl9yZXN1bHQkY2UsIGluZGV4LnJldHVybj1UUlVFKSRpeAptZWRpYW5faWR4IDwtIGFzLmludGVnZXIobGVuZ3RoKHNvcnRlZF9pZHhzKSAvIDIpICsgMQptYXhfaWR4IDwtIGFzLmludGVnZXIobGVuZ3RoKHNvcnRlZF9pZHhzKSkgLSAxCgpzZXFfbWVkaWFuX3ByZWRzIDwtIHN1YnNldChzZXFfcHJlZGljdGlvbnMsIHNlcV9udW0gPT0gbWF4X2lkeCkKCmJ4dCA8LSB1bmxpc3Qoc2VxX21lZGlhbl9wcmVkc1siWF9wcmVkX21lYW4iXSkKYnhzZXQgPC0gdW5saXN0KHNlcV9tZWRpYW5fcHJlZHNbIlhfcHJlZF92YXIiXSkKYnl0IDwtIHVubGlzdChzZXFfbWVkaWFuX3ByZWRzWyJZX3ByZWRfbWVhbiJdKQpieXNldCA8LSB1bmxpc3Qoc2VxX21lZGlhbl9wcmVkc1siWV9wcmVkX3ZhciJdKQoKTVJJbnB1dE9iamVjdCA8LSBtcl9pbnB1dCgKICBieCA9IGJ4dCwKICBieHNlID0gYnhzZXQsCiAgYnkgPSBieXQsCiAgYnlzZSA9IGJ5c2V0CikKTWVkaWFuRWdnZXJPYmogPC0gbXJfZWdnZXIoTVJJbnB1dE9iamVjdCkKYGBgCgoKYGBge3J9Cm1yX3Bsb3QoTVJJbnB1dE9iamVjdCwgbGluZT0iZWdnZXIiKQpgYGAKCgpgYGB7cn0KCmlmIChwYXJhbXMkc2F2ZV9wbG90cykgewogIHBuZyhmaWxlLnBhdGgocGFyYW1zJG91dHB1dF9kaXIsIG91dHB1dF9mbmFtZSksIHdpZHRoID0gNTEyLCBoZWlnaHQgPSAzMjApCn0KcGFyKG1mcm93PWMoMiwgMikpCgpmb3IgKHNlcSBpbiAoMTpuX3NlcXMpKSB7CiAgc2VxX2lfcHJlZGljdGlvbnMgPC0gc3Vic2V0KHNlcV9wcmVkaWN0aW9ucywgc2VxX251bSA9PSBzZXEpCgogIGJ4c2V0IDwtIHVubGlzdChzZXFfaV9wcmVkaWN0aW9uc1siWF9wcmVkX3ZhciJdKQogIGJ5c2V0IDwtIHVubGlzdChzZXFfaV9wcmVkaWN0aW9uc1siWV9wcmVkX3ZhciJdKQogIAogIHBsb3QoYnhzZXQsIGJ5c2V0LCBtYWluID0gc3RyX2MoIkNBIHZzLiBURiBzdGQgZXJyb3IiKSkKfQoKaWYgKHBhcmFtcyRzYXZlX3Bsb3RzKSB7CiAgZGV2Lm9mZigpCn0KYGBgCgpgYGB7cn0KZG1vZGUgPC0gZnVuY3Rpb24oeCkgewogIGRlbiA8LSBkZW5zaXR5KHgsIGtlcm5lbCA9IGMoImdhdXNzaWFuIikpCiAgKGRlbiR4W2RlbiR5ID09IG1heChkZW4keSldKQp9CmRtb2RlKGVnZ2VyX3Jlc3VsdCRjZSkKbWVkaWFuKGVnZ2VyX3Jlc3VsdCRjZSkKbGVuZ3RoKGVnZ2VyX3Jlc3VsdCRjZVtlZ2dlcl9yZXN1bHQkY2UgPCAwXSkKbGVuZ3RoKGVnZ2VyX3Jlc3VsdCRjZSkKYGBgCgpMZXQncyBhbHNvIGxvb2sgYXQgdGhlIGRpc3RyaWJ1dGlvbiBvZiBwbGVpb3Ryb3B5IHAtdmFsdWVzIGZvciBkaWZmZXJlbnQgc2VxdWVuY2VzLiBJZiB0aGV5J3JlIGNvbmNlbnRyYXRlZCBhcm91bmQgLjUsIHRoaXMgbWVhbnMgdGhhdCB0aGUgYWxnb3JpdGhtIGJlbGlldmVzIHRoZXJlIHRvIGJlIHZlcnkgbGl0dGxlIHBsZWlvdHJvcHkuCgpgYGB7cn0KcXBsb3QocGxlaW8sIGRhdGEgPSBlZ2dlcl9yZXN1bHQsIGJpbndpZHRoID0gLjA0KQpgYGAKCmBgYHtyfQpwIDwtIGdncGxvdDI6OmdncGxvdChlZ2dlcl9yZXN1bHQsIGFlcyh4ID0gY2UsIHkgPSBzdGQpKSArCiAgZ2VvbV9wb2ludCgpCnAKYGBgCgpgYGB7cn0KaWYgKHBhcmFtcyRzYXZlX3Bsb3RzKSB7CiAgb3V0cHV0X2ZuYW1lIDwtIHN0cl9jKGNlbGxfdHlwZSwgdGYsICJleF9lZmZfc2l6ZXMucG5nIiwgc2VwID0gIl8iKQogIG91dHB1dF9mcGF0aCA8LSBmaWxlLnBhdGgocGFyYW1zJG91dHB1dF9kaXIsIG91dHB1dF9mbmFtZSkKICBwbmcob3V0cHV0X2ZwYXRoLCB3aWR0aCA9IDUxMiwgaGVpZ2h0ID0gMzIwKQogIHBhcihtYXIgPSBjKDUsIDYsIDUsIDEpICsgLjEpCn0KCnBsb3QoYnh0LCBieXQsIHhsYWI9c3RyX2MoZXhwb3N1cmVfdGYsICIgRWZmZWN0IFNpemVzIiksIHlsYWI9c3RyX2Mob3V0Y29tZV90ZiwgIiBFZmZlY3QgU2l6ZXMiKSkKCmlmIChwYXJhbXMkc2F2ZV9wbG90cykgewogIGRldi5vZmYoKQp9CmBgYAo=